Epidemic modeling - Part 2

• Motivation for write-up
• Recall SEIR model equations
• Numerical solution to this deteministic population level model
  • Coding the SEIR model
  • Trying with average COVID-19 parameters:
• Effects of the different parameters on the SEIR simulation
  • Effect of $\beta$ ($r*\rho$)
  • Effect of $\sigma$ (average time from E → I)
  • Effect of $\gamma$ (average time from I → R)
• Discussion

Motivation for write-up

This is the 2nd part of a multi-part series blog post on modeling in epidemiology.

The COVID-19 pandemic has brought a lot of attention to study of epidemiology and more specifically to the various mathematical models that are used to inform public health policies. Everyone has been trying to understand the growth or slowing of new cases and trying to predict the necessary sanitary resources. This blog post attempts to explain the foundations for some of the most used models and enlighten the reader on two key points.

After introducing the concepts of compartmentalization and disease dynamics in the first blog post, this second part is focused on developing a deterministic numerical solution for the SEIR model discussed there.

Recall SEIR model equations

See the first blog post for derivation.

• Continuous-time:

  • $\frac{dS}{dt}=-\beta\frac{I(t)}{N}$
  • $\frac{dE}{dt}=\beta\frac{I(t)}{N} - \sigma E(t)$
  • $\frac{dI}{dt}=\sigma E(t) - \gamma I(t)$
  • $\frac{dR}{dt}=\gamma I(t)$

• Discrete-time:

  • $\Delta S = -\beta \frac{I}{N} * \Delta T$
  • $\Delta E = (\beta \frac{I}{N}-\sigma E) \Delta T$
  • $\Delta I = (\sigma E - \gamma I) \Delta T$
  • $\Delta R = \gamma I \Delta T$

Numerical solution to this deteministic population level model

Coding the SEIR model

#collapse_hide
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go

# Let's build a numerical solution
def seir_model(init, parms, days):
    S_0, E_0, I_0, R_0 = init
    Epd, Ipd, Rpd = [0], [0], [0]
    S, E, I, R = [S_0], [E_0], [I_0], [R_0]
    dt=0.1
    t = np.linspace(0,days,int(days/dt))
    sigma, beta, gam = parms
    for _ in t[1:]:
        next_S = S[-1] - beta*S[-1]*I[-1]*dt
        Epd.append(beta*S[-1]*I[-1]*dt)
        next_E = E[-1] + (beta*S[-1]*I[-1] - sigma*E[-1])*dt
        Ipd.append(sigma*E[-1]*dt)
        next_I = I[-1] + (sigma*E[-1] - gam*I[-1])*dt
        Rpd.append(gam*I[-1]*dt)
        next_R = R[-1] + (gam*I[-1])*dt
        S.append(next_S)
        E.append(next_E)
        I.append(next_I)
        R.append(next_R)
    return np.stack([S, E, I, R, Epd, Ipd, Rpd]).T

Trying with average COVID-19 parameters:

Parameters used for plot below:

• Days = 100
• Population = 10000
• Number of susceptible people on day 0 = 9999
• Number of exposed people on day 0 = 1
• No infected or recovered people on day 0
• $\sigma = 0.2$ (average of 5 days to go from exposed to infectious)
• $\beta = 1.75$ (average of $r=7$ contacts per day and $\rho = 25\%$ chance of a contact with an infectious person resulting in infection)
• $\gamma = 0.1$ (average of 10 days to go from infectious to recovered)

#collapse_hide
# Define parameters
days = 100
N = 10000
init = 1 - 1/N, 1/N, 0, 0
sigma = 1/5   # 1/5 --> 5 days on average to go from E --> I
beta = 1.75
gam = 1/10     # 1/11 --> 11 days on average to go from I --> R
parms = sigma, beta, gam

# Run simulation
results_avg = seir_model(init, parms, days)

#collapse_hide
fig = go.Figure(data=[       
    go.Scatter(name='S', x=np.linspace(0,days,days*10), y=results_avg.T[0], line={'dash':'solid', 'color':'blue'}),
    go.Scatter(name='E', x=np.linspace(0,days,days*10), y=results_avg.T[1], line={'dash':'solid', 'color':'yellow'}),
    go.Scatter(name='I', x=np.linspace(0,days,days*10), y=results_avg.T[2], line={'dash':'solid', 'color':'red'}), 
    go.Scatter(name='R', x=np.linspace(0,days,days*10), y=results_avg.T[3], line={'dash':'solid', 'color':'green'}),
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
  title={
      'text':'Deterministic SEIR model - COVID-19 parameters',
      'x':0.5,
      'xanchor':'center'
    }
)

fig.show()

Effects of the different parameters on the SEIR simulation

Effect of $\beta$ ($r*\rho$)

Let's have a look at the effect of $\beta$ on the SEIR simulation.

A higher $\beta$ can either mean a higher average number of contacts per day ($r$) in the population and/or a higher probabilit of transmission of disease from I → S.

The opposite holds also.

#collapse_hide
## Let's try to see how the model changes 
days = 180
N = 10000
init = 1 - 1/N, 1/N, 0, 0
sigma_avg = 1/5
beta_avg = 1.75
beta_noepi = 0.1
beta_low = 0.3
beta_high = 4
gam = 0.1     # 1/10 --> 10 days avergae I --> R
parms_avg = sigma_avg, beta_avg, gam
parms_noepi = sigma_avg, beta_noepi, gam 
parms_low = sigma_avg, beta_low, gam
parms_high = sigma_avg, beta_high, gam

# Run simulation
results_avg = seir_model(init, parms_avg, days)
results_noepi = seir_model(init, parms_noepi, days)
results_low = seir_model(init, parms_low, days)
results_high = seir_model(init, parms_high, days)

#collapse_hide
fig = go.Figure(data=[    
    go.Scatter(name=r'$S:\beta_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[0], line={'dash':'solid', 'color':'blue'}, legendgroup="COVID"),
    go.Scatter(name=r'$E:\beta_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[1], line={'dash':'solid', 'color':'yellow'}, legendgroup="COVID"), 
    go.Scatter(name=r'$I:\beta_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[2], line={'dash':'solid', 'color':'red'}, legendgroup="COVID"),
    go.Scatter(name=r'$R:\beta_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[3], line={'dash':'solid', 'color':'green'}, legendgroup="COVID"),
    go.Scatter(name=r'$S:\beta_{noepi}$', x=np.linspace(0,days,days*10), y=results_noepi.T[0], line={'dash':'dashdot','color':'blue'}, legendgroup="noepi"),
    go.Scatter(name=r'$E:\beta_{noepi}$', x=np.linspace(0,days,days*10), y=results_noepi.T[1], line={'dash':'dashdot', 'color':'yellow'}, legendgroup="noepi"),
    go.Scatter(name=r'$I:\beta_{noepi}$', x=np.linspace(0,days,days*10), y=results_noepi.T[2], line={'dash':'dashdot', 'color':'red'}, legendgroup="noepi"),
    go.Scatter(name=r'$R:\beta_{noepi}$', x=np.linspace(0,days,days*10), y=results_noepi.T[3], line={'dash':'dashdot', 'color':'green'}, legendgroup="noepi"),
    go.Scatter(name=r'$S:\beta_{low}$', x=np.linspace(0,days,days*10), y=results_low.T[0], line={'dash':'dash','color':'blue'}, legendgroup="low"),
    go.Scatter(name=r'$E:\beta_{low}$', x=np.linspace(0,days,days*10), y=results_low.T[1], line={'dash':'dash', 'color':'yellow'}, legendgroup="low"),
    go.Scatter(name=r'$I:\beta_{low}$', x=np.linspace(0,days,days*10), y=results_low.T[2], line={'dash':'dash', 'color':'red'}, legendgroup="low"),
    go.Scatter(name=r'$R:\beta_{low}$', x=np.linspace(0,days,days*10), y=results_low.T[3], line={'dash':'dash', 'color':'green'}, legendgroup="low"),
    go.Scatter(name=r'$S:\beta_{high}$', x=np.linspace(0,days,days*10), y=results_high.T[0], line={'dash':'dot', 'color':'blue'}, legendgroup="high"),
    go.Scatter(name=r'$E:\beta_{high}$', x=np.linspace(0,days,days*10), y=results_high.T[1], line={'dash':'dot', 'color':'yellow'}, legendgroup="high"),
    go.Scatter(name=r'$I:\beta_{high}$', x=np.linspace(0,days,days*10), y=results_high.T[2], line={'dash':'dot', 'color':'red'}, legendgroup="high"),
    go.Scatter(name=r'$R:\beta_{high}$', x=np.linspace(0,days,days*10), y=results_high.T[3], line={'dash':'dot', 'color':'green'}, legendgroup="high"),
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
    title={
        'text':r'$\text{Effect of } \beta \ \text{on Deterministic SEIR model}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

We notice a few things from the plot above on the impact of $\beta$:

• The higher $\beta$ is:
  • the faster the epidemic seems to propogate in the population
  • the higher the peak of infected individuals seems to be (meaning a higher chance hospital resources will be saturated)
• $\beta$ also appears to affect the overall number of people infected over the course of the epidemic - at one point, and keeping all other parameters the same, we reach a point with $\beta < 0.1$ where no epidemic even occurs (more on this later)

Effect of $\sigma$ (average time from E → I)

Let's have a look at the effect of $\sigma$ on the SEIR simulation.

A higher $\sigma$ means shorter average time to go from E → I.

The opposite holds also.

#collapse_hide
## Let's try to see how the model changes 
days = 180
N = 10000
init = 1 - 1/N, 1/N, 0, 0
sigma_fast = 1   # 1 --> Average 1 day from E --> I (ressembles SIR model)
sigma_slow = 1/100 #10 days on average, twice as long as COVID-19
sigma_avg = 1/5
beta = 1.75
gam = 0.1     # 1/10 --> 10 days avergae I --> R
parms_fastEI = sigma_fast, beta, gam
parms_slowEI = sigma_slow, beta, gam
parms_avg = sigma_avg, beta, gam

# Run simulation
results_fastEtoI = seir_model(init, parms_fastEI, days)
results_slowEtoI = seir_model(init, parms_slowEI, days)
results_avg = seir_model(init, parms_avg, days)

#collapse_hide
fig = go.Figure(data=[    
    go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[0], line={'dash':'solid', 'color':'blue'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[1], line={'dash':'solid', 'color':'yellow'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[2], line={'dash':'solid', 'color':'red'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[3], line={'dash':'solid', 'color':'green'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[0], line={'dash':'dash','color':'blue'}, legendgroup="high", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[1], line={'dash':'dash', 'color':'yellow'}, legendgroup="high", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[2], line={'dash':'dash', 'color':'red'}, legendgroup="high", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastEtoI.T[3], line={'dash':'dash', 'color':'green'}, legendgroup="high", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[0], line={'dash':'dot', 'color':'blue'}, legendgroup="slow", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[1], line={'dash':'dot', 'color':'yellow'}, legendgroup="slow", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[2], line={'dash':'dot', 'color':'red'}, legendgroup="slow", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\sigma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowEtoI.T[3], line={'dash':'dot', 'color':'green'}, legendgroup="slow", hoverinfo='x+y'),
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
    title={
        'text':r'$\text{Effect of } \sigma \ \text{on Deterministic SEIR model}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

We notice a few things from the plot above on the impact of the average time from E → I:

• The shorter the time, on average, from E → I:
  • the faster the epidemic propogates in the population
  • the higher the peak of infected individuals will be (meaning a higher chance hospital resources will be saturated)
• However, the time from E → I has no impact on the total number of individuals infected over the entire time of the epidemic (all are enventually infected in any case with these parameters)

Effect of $\gamma$ (average time from I → R)

Let's have a look at the effect of $\gamma$ on the SEIR simulation.

A higher $\gamma$ means shorter average time to go from I → R.

The opposite holds also.

#collapse_hide
## Let's try to see how the model changes 
days = 300
N = 10000
init = 1 - 1/N, 1/N, 0, 0
sigma_avg = 1/5
beta = 1.75
gam_avg = 1/10       # 1/10 --> 10 days average I --> R
gam_slow = 1/100  # 1/30 --> 30 days average I --> R
gam_fast = 1.5     # 1 --> 1 day average I --> R
parms_fastIR = sigma_avg, beta, gam_fast
parms_slowIR = sigma_avg, beta, gam_slow
parms_avg = sigma_avg, beta, gam_avg

# Run simulation
results_fastItoR = seir_model(init, parms_fastIR, days)
results_slowItoR = seir_model(init, parms_slowIR, days)
results_avg = seir_model(init, parms_avg, days)

#collapse_hide
fig = go.Figure(data=[    
    go.Scatter(name=r'$S:\gamma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[0], line={'dash':'solid', 'color':'blue'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[1], line={'dash':'solid', 'color':'yellow'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[2], line={'dash':'solid', 'color':'red'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{COVID}$', x=np.linspace(0,days,days*10), y=results_avg.T[3], line={'dash':'solid', 'color':'green'}, legendgroup="COVID", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastItoR.T[0], line={'dash':'dash','color':'blue'}, legendgroup="fast", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastItoR.T[1], line={'dash':'dash', 'color':'yellow'}, legendgroup="fast", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastItoR.T[2], line={'dash':'dash', 'color':'red'}, legendgroup="fast", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{fast}$', x=np.linspace(0,days,days*10), y=results_fastItoR.T[3], line={'dash':'dash', 'color':'green'}, legendgroup="fast", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowItoR.T[0], line={'dash':'dot', 'color':'blue'}, legendgroup="slow", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowItoR.T[1], line={'dash':'dot', 'color':'yellow'}, legendgroup="slow", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowItoR.T[2], line={'dash':'dot', 'color':'red'}, legendgroup="slow", hoverinfo='x+y'),
    go.Scatter(name=r'$S:\gamma_{slow}$', x=np.linspace(0,days,days*10), y=results_slowItoR.T[3], line={'dash':'dot', 'color':'green'}, legendgroup="slow", hoverinfo='x+y'),
])

fig.update_layout(
    xaxis_title = 'Day',
    yaxis_title = 'Proportion of population',
    title={
        'text':r'$\text{Effect of } \gamma \ \text{on Deterministic SEIR model}$',
        'x':0.5,
        'xanchor':'center'
    }
)

fig.show()

We notice a few things from the plot above on the impact of the average time from I → R:

• The longer the time, on average, from I → R:
  • the faster the epidemic propogates in the population
  • the higher the peak of infected individuals will be (meaning a higher chance hospital resources will be saturated)
• As opposed to the time from E → I but similarly as $\beta$, the time from I → R has an impact on the total number of individuals infected over the entire time of the epidemic (with no epidemic if $\gamma > 1.75$ and all other parameters are kept constant)

Discussion

So we can see the time periods for E → I and I → R, along with the value of $\beta$ are critical components in how the model will react.

Notably, no epidemic occurs if $\beta < \gamma$ (or if $\frac{\beta}{\gamma} < 1$)

In fact, $\frac{\beta}{\gamma}$ has major implications for the model and is called the basic reproduction number $R_0$.

Hence, if $R_0 < 1$ no epidemic occurs, while the opposite implies an epidemic.

There are major flaws with this model however. While this model is deterministic and uses average time to model $\sigma$ and $\gamma$, this is a major flaw and does not represent the reality for most diseases.

Part 3 of this blog series will discuss this further.